Coevolutionary dynamics on scale-free networks 
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We investigate Bak-Sneppen coevolution models on scale-free networks with various degree expo- 
nents 7 including random networks. For 7 > 3, the critical fitness value f c approaches to a nonzero 
finite value in the limit N — + 00, whereas f c approaches to zero as 2 < 7 < 3. These results are 
explained by showing analytically f c (N) ~ Aj < (k + l) 2 >jv on the networks with size N. The 
avalanche size distribution P(s) shows the normal power-law behavior for 7 > 3. In contrast, P(s) 
for 2 < 7 < 3 has two power-law regimes. One is a short regime for small s with a large exponent 
t\ and the other is a long regime for large s with a small exponent T2 (n > r%). The origin of the 
two power-regimes is explained by the dynamics on an artificially-made star-linked network. 

PACS numbers: 87.10.+e, 05.40.-a, 87.23.-n 
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Bak and Sneppen (BS) Q has introduced an excellent 
model to explain the evolution of bio-species which ex- 
hibits the punctuated equilibrium behavior |2j. BS model 
has two important features, coevolution of the interacting 
species and the intermittent bursts of activity separating 
relatively long periods of the stasis. In BS model the 
ecosystem evolves into a self-organized criticality with 
avalanches of mutations occurring all scales. Aside from 
its importance for the evolution BS model has been also 
shown to have rich scaling behaviors |3j]. 

Since BS model was suggested, the model has been 
extensively studied on regular lattices or networks Q. 
However, many important bio-systems have been eluci- 
dated to form nontrivial networks by the recently de- 
veloped network theories Q. Important examples are 
metabolic network, cellular network, and protein network 
HL EJ El ■ Especially the important bio- networks are 
scale- free networks (SFNs) 0], in which the degree dis- 
tribution p(k) satisfies a power lawp(fc) ~ fc~ 7 [j]. Thus 
it is important to study the BS dynamics on SFNs or to 
find out how the base structure of interacting biological 
elements (cells, proteins, or species) affects the evolution- 
ary change or dynamics of the given bio-system. Until 
now BS models on the nontrivial networks were not inves- 
tigated extensively. Christensen et al. [9| have studied 
BS model on random networks (RNs). Kulkani et al. 
[Tl| studied BS model on small- world networks. Slania 
and Kotrla ^1 studied the forward avalanches of a sort 
of extremal dynamics with evolving networks. Moreno 
and Vazquez [lj| studied BS model only on a SFN with 
7 = 3. 

In this letter, we will study BS models on SFNs in 
complete and comprehensive ways. One of the main pur- 
poses of this study is to find which structure of interact- 
ing species is the most stable network or most close to 
mutation-free network under the coevolationary change 
with interacting species. As is well-known, SFNs with the 
degree exponent 2 < 7 < 3 are physically much different 
from those with 7 > 3 Q . We study BS models not only 
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on SFNs with 2 < 7 < 3 but also on SFNs with 7 > 3 
including random networks (or SFN with 7 = 00). As we 
shall see, two important results are found in this study. 
First, the critical fitness value f c of BS models for 7 < 3 
is shown to have the limiting behavior f c (N) — ► when 
the number of nodes N of the network goes to infinity. In 
contrast, f c approaches finite nonzero value as N — > 00 
for 7 > 3. Furthermore, f c (N) on SFNs with finite N is 
shown to satisfy the relation f c (N) ~ < ^°"^ >jv , which 
is also directly supported by simulation. Second, for 
2 < 7 < 3 the distribution of avalanches is shown to 
have two power-law regimes. To find the origin of this 
anomalous behavior of avalanches we also study BS mod- 
els on an artificially-made star-linked network and find 
the similar two power-law regimes. 

We now explain the model treated in this letter. All 
the models are denned on a graph Gr — {N,K}, where 
N is the number of nodes and K is the number of de- 
grees with the average degree < k >= 2K/N. Initially, a 
random fitness value € [0, 1] is assigned to each node 
i = 1, N. At each time step, the system is updated by 
the following two rules: (I) first assign new fitness value 
to the node with the smallest fitness value f m in- (II) 
Second assign new fitness values to the nodes which are 
directly connected to the node with f m in- We use SFNs 
with the various degree exponent 7 as Gr = {N, K}. To 
generate SFNs, we use the static model [lj] instead of 
preferential attachment algorithm Q. 

To understand the dependence of the critical fitness 
value f c (N) on 7, we generate SFNs with 7 = 00, 5.7 ~ 
2.15. To exclude the effects of finite percolation clusters 
9| and to see the effect of network structure itself, all the 
networks are made to have average degree < k >= 4. To 
understand the dependence on number of nodes N, the 
networks with the sizes N = 10 3 ~ I0 6 are generated for 
each 7. To determine the critical fitness value f c (N), we 
consider f m i n as a function of the total number of updates 
s 0. Initially, f m in(s = 0) is the gap G(0), where G(s) 
is the maximum of all fmin(s') for < s' < s 3]. When 
G(s) jumps to a new higher value, there are no nodes in 
the system with fi(s) < G(s). Thus lim^oo Gn(s) = 
UN). 
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FIG. 1: Semi-log plot of the threshold fc(N) versus 1/N on 
RN and on SFNs with 7 = 5.7, 4.3, and 3.5. Used networks 
sizes for each network are TV = 10 3 , 10 4 , 10 5 , and 10 6 . The 
solid lines between data points are obtained by simple linear 
interpolations. 



We measure f c {N) on the various SFNs. Fig. 1 shows 
the plot of f c (N) versus 1/N for SFNs with various 7. 
The values of critical fitness f c {N — > 00) evaluated from 
data in Fig. 1 are 0.21(1), 0.19(1), 0.15(1), and 0.09(1) 
for 7 = 00, 5.7, 4.3, and 3.5. The results in Fig. 1 mean 
that for 7 > 3, f c {N — * 00) — ► const. (> 0). 

Fig. 2 shows the plot of f c (N) versus 1/N for 2 < 
7 < 3. For 7 = 3, fc(N) nicely satisfies the relation, 
f c (N) ~ 1/lnTV [13. For 2 < 7 < 3, / c (iV)'s seem to 
follow a power-law f c (N) ~ N~ v and approach to zero 
as N goes to 00. In contrast to the results in Fig. 1, 
f c -> for 2 < 7 < 3. 

In the RN, every pair of nodes are randomly con- 
nected and the degree distribution is a Poisson distri- 
bution HI!. So the BS model on RN is a good re- 
alization of the mean-field-type random neighbor model. 
In the random neighbor model, the fitness values of the 
randomly selected (m — 1) nodes as well as the node 
with f min are updated and f c = 1/m [Tjj. The result 
/ c (oo) = 0.21(1) on RN is very close to 
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which is expected one from the random neighbor model 
by setting <J;>+l = m|. In the steady state of BS 
model, the probability measure P(f < f c ) is 0. Suppose 
the case that the number of updates for each step is fixed 
as to, as in the random neighbor model. To sustain the 
steady state in the case, at most one new fitness value 
should be less than f c and the other m — 1 new values 
should be larger than f c ^3 • Therefore we can easily see 
mf c = 1 or f c = 1/m. 

On a network the number of updats depends on the 
degree of the node with f m i n and the probability which 
a node with degree k is connected to the node with f m in 
should be proportional to k. For an updating step the 
probability that a node with degree k is updated is pro- 
portional to k + 1, because the node itself can be the 



node with / TO , n . Therefore, after an arbitrary update, 
the probability P m i n {k) of a node with degree k being 
the node with f m in is proportional to k + 1 . This means 
that P m in{k) in the steady state should be proportional 

to k + 1, or P mm (k) = ^ffffiw = <^+T (fc + 1)p(fc) - 
The average number N up date of the nodes updated for 
one updating process is therefore 
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When the number of updates is fixed as to, Eq. |j5J re- 
produces the mean-field result f c = 1/m. In SFNs with 
p(k) ~ fc~ 7 , Eq. @ becomes 
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2 < 7 < 3. 
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Eq. J2J explains the results in Figs. 1 and 2 including 
the result f c ~ for 7 = 3. For 2 < 7 < 3, measured 
f c (N) is fitted to the relation f c (N) = A/ < k 2 > N , 
where A is constant and < k 2 >n is < k 2 > for the 
network with the size N. The fitted lines in Fig. 2 show 
that the relation f c (N) = A/ < k 2 >n holds well and 
directly supports Eq. @. 

An avalanche in Bak-Sneppen model is defined as the 
sequential step s for which the minimal site has a fitness 
value smaller than given f Q Q. For each network, we 
choose f to satisfy (f c (N) - fo)/fc(N) = 0.05. The 
probability distribution P(s) of avalanche size s on the 
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FIG. 2: Log-log plot of f c (N) and A/ < k(N) 2 > N versus 
1/N on SFNs with 7 = 2.75, 2.40, and 2.15. Symbols are 
for fc(N) and the lines are for A/ < k(N) 2 >jv, where A 
is a constant. The top inset shows the plot of f c (N) versus 
1/lnN for 7 = 3.0. 
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TABLE I: Two power-law exponents, n and T2 for SFNs with 
7 < 3. 
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FIG. 3: Log-log plot of the avalanche size distribution P(s) 
on SFNs with 7 = 5.7, 7 = 4.3, 7 = 3.5 and on RN (Inset). 
The curves for 7 = 5.7, 7 = 4.3 and RN denote the fits 
of the form P(s) — As~ T exp( — s/s c ) to the data. Obtained 
exponents are r = 1.5 for both 7 = 5.7 and RN, and r = 1.65 
for 7 = 4.5. The line for 7 = 3.5 denotes the fit of the form 
P(s) = As~ t {t = 1.65) without cutoff. 



networks with the size N = 10 6 are shown in Fig. 3 and 
Fig. 4. All the data in Figs. 3 and 4 are taken in the 
steady-states. 

As is shown in Fig. 3, P(s) in SFNs with 7 > 3 in- 
cluding RN satisfy the normal power-law behavior with 
an exponential cutoff as P(s) = As~ T cxp(— s/s c ). The 
curves in Fig. 3 represent the fitted curves to data for 
P(s). From those fittings the obtained values for r are 
1.5 for RN and 7 = 5.7, and 1.65 for 7 = 4.3. The re- 
sult for RN and SFN with 7 = 5.7 is expected from the 
random neighbor model [l(|. As 7 decreases to 4.0 or 
so t increases to 1.65. For 7 = 3.5, however, the best 
fitting function is P(s) = Bs~ T with r = 1.65 and we 
cannot find the cut-off-dependent behavior within our 
data. Instead, it is even observed that tails of measured 
data for 7 = 3.5 around s — 10 3 seem to deviate from 
the fitting function P(s) = Bs~ T and are lager than val- 
ues estimated from the best fitting function. This rather 
anomalous tail behavior of P(s) for 7 = 3.5 should be the 
signal of the anomalous behavior of P(s) for 2 < 7 < 3. 

In contrast to the simple power-law behavior for 7 > 3, 
anomalous behavior for P(s) shows up for 2 < 7 < 3 
(Fig. 4). We can see two power- law regimes clearly for 
P(s) in Fig. 4. Initially the avalanche size distribution 
follows P(s) ~ s~ Tl about 1 decade or so. After this 
short initial power-law regime, the long second power- 
law regime appears as P(s) ~ s _T2 , where t\ > T2- The 
measured exponents 77, T2 are summarized in Table I. 

Compared to the behavior of the avalanche size dis- 
tribution for 7 > 3, this anomalous behavior of P(s) is 
very peculiar. In the steady state, it is expected that the 
node with / TO , n (the minimal node) is most frequently 
found among the last updated nodes ^3 an d then the 
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minimal node locally performs a random walk. However, 
there can be longer jumps of any length with a very low 
probability. If this kind of a jumpy random walk is the 
motion of the minimal node, then a subnetwork consists 
of a hub node (center node) and many slave nodes di- 
rectly linked to the hub should be important to decide 
the behavior of P(s). Due to the jumpy random walk 
behavior, the more slave nodes the hub node has, the 
longer stay of the minimal node or the longer avalanche 
exists at the given subnetwork. This effect explains the 
second power-law regime with the exponent t-i in Fig. 
4, because < k 2 > diverges for 2 < 7 < 3, and so the 
subnetwork of a hub node and many slave nodes should 
be the main substructure in SFNs with 2 < 7 < 3. Evi- 
dently, the jumpy steps of the jumpy random walk make 
the shorter avalanches possible and this effect explains 
the first power-law regime with the exponent 17 . 

To support the qualitative explanation of the two 
power-law regimes, we consider an artificially-made star- 
linked network shown in Fig. 5. In the star-linked net- 
work, a main subnetwork consists of a center (star) node 
and many dangling slave nodes linked directly to the star 
node. Then the center nodes are linked hierarchically to 
one after another as sketched in Fig. 5(a). We make 
a star-linked network in which there are 25 base sub- 




FIG. 4: Log-log plot of P(s) on SFNs with 7 = 3 (top inset), 
2.75, 2.4, and 2.15. Two crossing lines for each data sets 
denote the two power-law regimes, P(s) = As~ T1 and P(S) = 
Bs~ T2 . Obtained exponents, ri and T2, are shown in Table I. 
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FIG. 5: (a) Schematic diagram of a star-linked network which 
consists of 25 subnetworks with 500, 480, and 20 dangling 
slave node, (b) Plot of P(s) on the star-linked network struc- 
ture. Two power-law regimes with P(s) = As~ T1 (n = 3.7) 
and P(s) = Bs~ T2 (t2 = 1.27) are clearly shown by the lines. 

networks with 500, 480, and 20 slave nodes, respec- 
tively. In this network, we perform BS dynamics and 
find f c = 0.123. P(s) is also measured on the star-linked 
network and is shown in Fig. 5(b). We find the very 
two power-law regimes with the exponents t\ = 3.7 and 
T2 = 1.27. The plateau between two power-regimes in the 
data of P(s) in Fig. 5(b) is probably from the discrete 
distribution of the number of slave nodes. 

In conclusion, we study BS models on SFNs with var- 
ious 7. For 7 > 3, f c approaches to a nonzero value 
in the limit N — > 00 and P(s) shows normal power-law 
behavior with r > 1.5. For 7 < 3, f c approaches to 
zero as f c (N) ~ A/ < K 2 >n and P(s) has two power- 
law regimes. The origin of the two power-regimes are 
explained by the dynamics on a star-linked network. 

In Ref. 13], BS dynamics only on a SFN with 7 = 3 
was studied and the only meaningful numerical result 
was to show f c (N) ~ 1/lniV. Ref. ,13] suggested a re- 
lation similar to Eq. from a rate equation which was 
obtained by a naive and immature analogy of BS dy- 



namics to the epidemic dynamics on SFNs 15]. However 
the rate equation should never be the exact one. Even 
the exact rate equation for the simple random neighbor 
model [13 is much more complex than that of Ref. 0] 
or the epidemic dynamics. The correct rate equation 
for BS dynamics on SFNs must be derived by consider- 
ing all the terms of the rate equation in Ref. ^(j an d 
the base network structure simultaneously and correctly. 
The derivation of the correct rate equation should be 
a subject for the future study. In Ref. 0] they argued 
P(s) for 7 = 3 satisfies a simple power-law with r ~ 1.55. 
By the brute-forced fit of the relation P(s) ~ s~ T to our 
data in Fig. 4, we also obtain r ~ 1.6 for 7 = 3. However, 
this blind application of the simple power law should be 
wrong and there should exist the two-power law regimes 
even for 7 = 3. One can easily identify the two power-law 
regimes in the P(s) data of Ref. [jjj rather clearly al- 
though the tail parts of their data are qualitatively poor 
and show large fluctuations. 

The occurrence of two power-law regimes for P(s) was 
also found in BS dynamics on small- world networks 111! 
and in an extremal dynamics with evolving networks [12j . 
However the origins of the two power-law regimes were 
completely different from ours. The origin in the small- 
world networks was argued to be the long range con- 
nectivity of the networks [Tl| . The extremal dynamics 
with evolving random networks |12| changes the network 
structure and is not exactly the same as BS dynamics. 
Furthermore the evolving network develop many discon- 
nected clusters. In the model the forward avalanches 
are mainly measured. The forward avalanches 12] should 
be affected by the dynamical aggregate and splitting of 
subnetworks by the extremal dynamics, which should be 
the origin of the two power-law regimes. In contrast 
our avalanches of BS dynamics is measured on a fully- 
connected static scale-free network and should not be di- 
rectly comparable to the avalanches on dynamically vary- 
ing networks. 
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